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SECTION I 
INTRODUCTION 


Ever since the earliest days of manned flight, panel flutter has been 
known as one of the most important problems in the design of aircraft, 
missiles. launched vehicles, and spacecraft. Extensive progress in wind 
tunnel tests and theoretical studies has been achieved. The basic theories 
and an account of the developments on panel flutter can be found in, amona 
other books, a recent text by Dowell (Reference 1). / list of keyed 
bibliography and collection of some significant survey papers and original 


papers were prepared by Garrick (Reference 2). 


The theoretical solution for a panel flutter problem usually requires 
an accurate aerodynamic and structural theory to formulate a set of complex 
eigenvalue equations of motion interacted between the panel and the flow. 
One of the most common theoretical methods is the modal method where the 
aerodynamic pressure and the inertial and elastic forces of the panei are 
obtained by assuming the displacements as composed of a number of natural 
modes and generalized coordinates. The natural frecuenctes and corresponding 
normal mode shapes are obtained either theoretically or experimentally, 
or both. Since the finite element method is powerful and practical in the 
free vibration analysis of panels with arbitrary geometrical anc boundary 


conditions, it is commonly used in obtaining natural frecuencies and node 


shapes in the modal method. 


fc an alternative approach, the finite element workers have 


formulated the matrix of aerodynamic pressure by using the 


AT ma te ta 


displacement functions as composed uf the nodal degrees of freedom and shape 
functions rather than the generalized coordinates and natural mode shapes. 
Such approach can directly solve for the flutter frequencies and corresponding 
normal mode shapes without having to seek the natural frequencies and modes 
and choose the number of modes before the eigensolution, and also without 
having to compute the flutter mode shapes after the eigensolution, Such 
approach permits expression of the equations of motion in an elegant and 
straightforward form. Such approach also permits generality in panel 
configurations and boundary conditions, and allows for flexibility to 


accurately include physical effect such as in-piane forces. 


The finite element method was first extended to the panel flutter 
problems by Olson (Reference 3). He formulated the aerodynamic matrix 
explicitly for an infinite plate element. Olson (Reference 4) later 
formulated the aerodynamic matrices for two rectangular (12 and 16 d.o.f.) 
and an 18 d.o.f. triangular plate elements. Simultaneously, but 
independently, Appa and Somashekar (Reference 5) formulated the aerodynamic 
matrix for a 12 d.o.f. rectangular plate element. Appa, Somashekar and 
Shah (Reference 6) later extended their work by accounting for skew panels 
and yawed flow by means of coordinate transformation. Sander, Bon, and 


Geradin (Reference 7) employed the CQ conforming quadrilateral plate finite 


element for fluttcr analysis of rectangular panels with yawed flow and 


in-plane stresses. 


In all the above finite element works (References 3 - 7), Lighthill's 
linearized piston theory was employed. The Mach numbers considered were 
above approximately 1.6. Recently, Yang (Reference &) developed a finite 
element procedure using the exact linearized two-dimensional theory (strip 
theory) for unsteady supersonic flow to formulate an infinite plate finite 
element by means of numerical integration. The flutter speed considered 
thus could be in the lower supersonic range. Such formulation cannot, 
however, be adequately applied to the more general case of rectangular 


panels with finite aspect ratios. 


In this report, the three-dimensional supersonic unsteady potential 
fiow theory ts employed to formulate the plate finite elements so that the 
flutter problems of the finite panels in low supersonic range can be 


treated. 


The aerodynamic matrix is derived by using the principle of virtual 
work. The aerodynamic forces or velocity potentials that produce the work 
are obtained for each d.o.f. by the Mach box method. fach finite element 
is divided into several boxes, The aerodynamic influence coefficients 
for each pair of sending and receiving boxes are evaluated, for cach d.o.f., 
by the method cf Gaussian quadrature. The velocity potential at each 
receiving box is obtained, for cach d.o.f., by summation of the product 
of corresponding downwash and influence coefficients for all sending boxes. 
]t should be noted that a similar box method was used by Cunningham (Reference 9) 
in conjunction with a Galertin's mudal wethod for panct flutter analysis, 
Such 3-D aerodynamic theory was also employed by Dowell and Voss (Reference 10) 


in a theoretical and experimental correlation study of panel flutter. 


bo 


The 16 d.o.f. conforming rectangular plate element (Reference 11) 
was used for example demonstrations. Flutter boundaries were found for 
clamped rectangular panels with various aspect ratios. The thickness 
ratios required to prevent flutter of an aluminum panel at sea level were 
plotted for Mach numbers ranging from 1.05 to 3. Results found are in 


favorable agreement with Cunninhgam's solution (Reference 9). 


The initial in-plane tensile stresses were then included in finding 
the flutter boundaries for a clamped aluminum panel. The thickness ratios 
required to prevent flutter of the panel at 25,000 feet altitude were 
obtained for Mach numbers ranging from 1.05 to 2.0 for various values of 
tension. It was found that the dominating flutter mode changed abruptly 


as Mach number was varied. 


SECTION I] 
FORMULATION OF EQUATIONS OF MOTION 


The free vibration equations of motion for a finite element panel 
subjected to the effect of stiffness, in-plane force, inertia, and 


aerodynamic pressure may be wricten as 


[K]{q} + [N}{q} + [M}{q) + [A]{q} = {0} (1) 


where [K], [N], [M], and [A] are, respectively, the stiffness, incremental 
stiffness, mass, and aerodynamic matrices assembled for the whole finite 
element system. The vector {q} contains the nodal degrees of freedom for 
the whole panel system. In this section, only the system aerodynamic 


matrix is formulated. 
1. PRINCIPLE OF VIRTUAL WORK 


For a system of plate finite elements, the deflection and aerodynamic 
pressure may be written by separating the time and space variables as, 
- iwt 
z(x,y.t) = z(x,yJe 


B(x,yst) = plx.y)e'%™ 


where the coordinates are defined in Fig. 1. 


The strain energy for the panel system is equal to the work produced 


by the aerodynamic pressure 


{| z(xsy)p(xsy) dxdy (3) 


Figure 1. Panel diviced into tinite clements, (solic Lines) 
and boxes (dash lines) with coordinates, Aimeansions, erd Mach cane 
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Figure 2. Three tones af integrals for cumpeat ing the 
aovexignamicG dp fluence coef iedents 


The deflection function for the total finite element system may be assumed as 


N 
z(x,y) = } f (y)q, = cf}! tq} (4) 
jr 


where f (xy) represents the assembled shape function corresponding to the 
nodal d.o.f. as and N is the number of degrees of freedom for the entire 
panel system. Usually, 95 represents deflection, slopes, twist, and curvatures 


of the plate at the nodal point j. 


Similarly, the aerodynamic pressure may be assumed as 
= y)q, = (Pig) (5) 
p(xsy) =) Py ly)ay + {q 


where PCy) is defined in Eq. (8) as the pressure function associated with 


the shape function f Gay) and d.o.f. "a5" 


Substituting Eqs. (4) and (5) into (3), the strain energy expression 


becomes 


Rh] — 


U= + igh [Al(a) (6) 


with 


yY dxdy (7) 


[A] {| FHP 


Using the principle of virtual work, it can be interpreted that Eq. (7) 


yields the aerodynaniic matrix, 
Z@. AERODYNANIC PMTRIX AND AERODYNAMIC PRESSURE 


The acrodynamic perturbation pressure is obtained from the linearized 


three-dimensional supersonic unsteady potential flow theory. When associated 


with the d.o.f. nq and the shape function fQGy), the perturbation 
pressure may be expressed in terms of the velocity potential ttirough the 
relation (see, for example, Ref. 12), 


36 - 
MI Gee a.) (8) 
& 3x 


P.(x,y) = 
5 y) V 


where 5 is the velocity potential associated with d.o.f. "q.". 


The coefficient for the ith row and jth column of the aerodynamic 


matrix [A] is obtained by substituting Eq. (8) into Eq. (7). 


i 
cd. 


ag 7 EYP Qs a a veomatuy) (9) 


ax 
Performing integration by parts and assuming restrained trailing edge, 


fq. (9) becomes 


of. 
ms cies: Oe al fay 7 (10 
Aj. = = {| ails 1 yr fi )dexddlwy) (10) 


where the quantity within the parentheses is the complex conjugate of the 
downwash ratio for unit d.o.f. "a," and unit panel length. The coefficient 


Aj is evaluated numerically through the following simple summation 


ON A GU ge hata ez 
A.. = = y SP i§ f.) (Box Area) (11) 


where the summation is to be mau for every box and the quantities related 
to tj anu fi are evaluated numerically at the center of each box. The method 
for numerically evaluating the velocity potential at the center of each box is 


given as follows. 


& 


3. VELOCITY POTENTIAL AND AERODYNAMIC INFLUENCE COEFFICIENT 


The method of Mach box as employed by Cunningham in Reference 9 is 
used here to derive the aerodynamic influence coefficients between each 


pair of boxes and the resulting velocity potentials at each box. 


In this method, each finite element is divided ints several equal-size bexes 
as shown in Fig. 1. The numbers of boxes in the stream and cross-stream directions 


are defined as BY and By ss respectively. The widtn and length of each box are 


defined as ce and te, resnectively, where e = w/B. and t = 2B gf ib . The boxes 


X 


G an oe - ee Saux h Ma oes Sratle. oe atta, Sa, ia! aes 2 a haw 
are tumbered in sequence fron the origin with ii (or 4) and a (or v) aS thé box 


~v 


number in the stream and the cross-stream directions, respectively. The numbers 
(m,n) and (+,v) refer to the receiving and sending boxes, respectively. 

The boxes are assumed as sufficiently small so that the downwash over any 
sending box is considered as uniformly distributed at any instant, and the 
resulting perturbation pressure et the center of each receiving box represents 
the average of the pressure distributed over the box. 

The velocity potential at the center of a receiving box (m,n) due to a 
uniformly distributed hut otherwise unspecified downwash wO>,.) over the sending 


box (A,v) can be expressed, for harmonic motion, as 


¢(m,n) = « widX,v)a, (r,s) (17) 


¢ 
where the relative locations in stream and cross-stream directions, respectively, 
between the two boxes are defined as r-7 ti - 4 ands tn - vy. The aerodynainic 

ry 


4 Auer eh Ru 8 Wh i eh nN fe 
Vey Uoree CULTFEICTOHE trum nererenbe 9 Ia, 


Bo 2 4 45 
Ube ep ase coe ay es Vv> cos[sr(u2-v")*] - 1 
a (r,s) = 3. | e W3E/ {cos i mika cos! ie | tt dvd 
® ae St Vi (ue - v°)? 


where all the parameters are defined in the List of Symbols. 
In Eq. (13), the surface integration limits U;, Us, Vy}, and Va result in 


only three different forms inp Loos and In, as shown in Fig. 2. 


The first funn is for any portion of a sending box(*, v) cut by both sides of 


7 

the Mach coné so that v» = -vi = u ands = 0. 
be LT? gTCe/7 tg. (isa (14) ! 
Gi 2 2M 


where Jg is the Bessel function of the first kind of order zero. 
The second form is for portions of a box cut by one side of the Mach cone so 


that the limits v, = u end v, = 2s - 12> 1. 


aa, 4 a\h 
We Pee aly ob costae SV ya) 1 
In, = i. f e 4 (85./ {cos 1 + | eee Ramo ————— dee (75) 
, ) "4 (u? - y?y2 J 
The third form is for boxes that are completely within the Mach cone and 
also for portions of boxes ahead of the point where the Mach line cuts the side 


1 of the box, 


ise ahr cello 
Oy ot oe ok ys Vz COstafueeve )tp = 
les * Z : Teel og 1 Tiss | 2H acevo avy 
, uy vy (ue - ve}? 
(16) 
{ 
The complete «,(r,s) for any one sending box (s,) consists of loys los or loys 


or a comuination cof t,. and L.., or of |, . and | 
Qa. G: hh. G 


The types of integrals and 


limits of inteyration for conputing «.(r,5) tor all possible relative Jocalions 


of box(4, vyand the Mach cone from(s, n)are given in Reference 9. 


10 


———— eee 


The evaluation of the above three integrals is carried out through the 
method of Gaussian quadrature. In the subsequent numerical examples, three 
Gaussian aout we used in both x and y directions for computing In 5: Five 
Gaussian points in both x and y directions are used for computing ley and Io: 

Once all possible values of the aerodynamic intluence coefficient a (78) 
are obtained for a certain shape function f OGY), the total velocity potential 


at the center of a receiving box (m,n) for the downwash associated with the sanie 


shape function is a weighted sum of the (m,n) defined in Eq. (72) 
(rss) (17) 


The summation is extended over all the sending boxes. The downwash ratios 
We Osv)/V for a unit d.o.f. "a," are the total time derivatives of the shape 
function f(x ye 
W. af. : 
peeps i gery) (18) 
ax 2 
Once the velocity potentials are obtained for each box for each shape function, 


they can readily be substituted into Eq. (11) for computing the aerodynamic 


matrix for the entire panel system. 


SECTION III 


FORMULATION FOR FLUTTER DETERMINANT 
IN TERMS OF FINITE ELEMENT SHAPE FUNCTIONS 


In the present finite element flutter formulation, the parameters are 
so grouped and nondimensionalized that the solution is general enough to 
include every parameter. Assuming harmonic motion with natural frequency w, 


the equations of motion (1) are rewritten in an element form as 


fo fe,k} + Festn}} - featm) - eyfA]]}ta) = (0) 


where % = wi7/w2(1 + ig) is the flutter eigenvalue parameter, and 


Cc 
g ~ 
cy = u(t/a)"d 7 C5. = otayar? : 


w)2ahe" 
2 ; 2 by3 € a 2 (tel 
c3 = u(2/a) Gy) k,? TE, 
The element matrix terms are obtained as 
2 a2" 2 2 2 
- 1 /-l off. a¢f ach h) f. o¢f. ay? ant act 
Ki 7 ax ge Sy Noes gee aad) ar ae 
J 70-0) ox ox oy ay” 9X ay 
a7f. a2f. > of. 37f, 
+ v(2)? 2 + 201 - v)(®)° + diay 
b =2 -2 b =.- 
ox 80ooy OXay dXdy 
1 1 oof, af. N of. of. N . af. of. N of. of. 
face | | [ae ee ge Ss 
VJ G40 ax ax x dy ay X 3% dy X 9x dy 
VA _ 
eo f.f. dxd 
Mj : ly 
>. 3 
PACs eee © cee erm ae 
AG Lye vty) 
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(19) 


(20) 


The shape functions fi (xs¥) are usually cubic or higher order functions of 
nondimensional coordinate variables x and y. Their differentiations with 
respect. to xX and y are performed analytically. The subscript "i" indicates 
the degree of freedom number with which the shape functien is associated. 
Eqs. (19-24) are suitable for any plate finite element so long as it is a 
displacement model with assumed shape functions. According to the state-of- 
the-art of the finite element development, the stiffness matrix [k], mass 
matrix [m], and incremental stiffness matrix [n] have been formulated 
analytically and explicitly for almost every plate and shell finite element. 
Only the aercdynamic matrix [A] remains to be formulated and it is to be 


obtained by numerical integration here. 


Eq. (19) constitutes an eigenvalue problem. The flutter solution is 
obtained by first assuming a valuc of reduced frequency Koy then varying 
the air-panel mass ratio 1/y incrementally and solving for the corresponding 
eigenvalues © or (w /0)70 + ig). When the structural damping coefficient 
g changes its value from negative to positive, the panel goes from the stable 
region to unstable region, and vice versa. The values of k, and 1/p that 
correspond to zero g value define the flutter boundary and the corresponding 


mode shape defines the flutter mode. The complex eigenvalue problem is solvec 


by uSing a subroutine in EISPACK provided by Argonne National Laboratory. 


Before performing analysis for each class of problems, a convergence study 
by varying the meshes of Mach boxes and finite elements must be made in order 
to find the suitable meshes needed for obtaining converged results. Such meshes 
depend on the panel geometry, boundary conditions, flow speed, and flutter modes. 
However, due to the enormous computations needed for obtaining massive data in 
this study, not the finest meshes are used. Most of the present results are 


compared with those obtained by the modal method (Reference 9) and good 


results are obtained. 
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SECTION IV 


RESULTS 


The 16 d.o.f. conforming rectangular plate finite elements were 
employed to demonstrate the present formulation and procedure. The examples 
chosen were rectangular panels, stressed as well as unstressed, with clamped 
edges. For the unstressed panels, a sophisticated analytical solution 
by Cunningham (Reference 9) using the box method (400 boxes) and Galerkin's 


modal approach (6 to 16 modes) was available for comparison. 


In all examples studied here, the flutter mode shapes were assumed as 
symmetrical about the center chord-line of the panel. Thus, only half of 
each panel was modeled by finite elements. In all cases, two elements in 
the cross stream direction were used for half of the panel. The number of 
elements in the stream direction varied from 4 to 10 dependent on the 
dominating flutter mode shapes. Each element was divided into 4 by 2 boxes 


in the stream and cross stream directions, respectively. 


It is important to note that symmetry does not exist for the mesh of 
boxes unless the Mach cone apex is located at the center chord-line of tie 
panel. The boxes on the other (disregarded) half of the panel can, however, 


still be accounted for since their deflection shapes are known by symmetry. 


Case 1 Clamped Panels with Various Chord-span Ratios at M = 1.3 


The rectangular panels with all edges clamped and chord-span ratios 


(¢/w) equal to 0, 1/4, 1/2, 1, 2, and 4 were studied for M = 1.3. 
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The resulting flutter boundaries are presented as plots of 
each dominating mode with the air-panel mass ratio 1/, as the vertical 
coordinate and the stiffness parameter wi tiN as the horizontal coordinate. The 
values of reduced frequency Ke or ws/V are marked along each curve. The 
result ‘c: the infinite panel with s/w = 0 are in exact agreement with those 


obtained by Cunningham, therefore, they are not presented here. 


Fig. 3 shows that for the wide panel with &/w = 1/4, the first mode 
flutter boundary crosses the second mode flutter boundary. The critical 
flutter boundaries are thus dominated by the first mode in the high mass ratio 
range and by the second mode in the low mass ratio range. A slight amount of 
structural damping was also included (g = 0.01). The results agree well 


with those obtained by Cunningham (Reference 9). 


The first flutter mode shape (real part) corresponding to point A (ky = 0.635) 
in Fig. 3 is shown in Fig. 4. It is seen that, due to the effect of the 
flow, the mode is not symmetrical about the cross stream centerline. The 


maximum deflection occurs at a point behind the center of the panel. The 


second flutter mode shape (real part) for the center chord-line of the panel 


corresponding to puint B (kK, = 1.59) in fig. 3 is shown in Fig. 5. The effect of the 


flow that pushes down the front part of the panel is seen. 


Fig. @ shows that, as the panel width is reduced (*/w = 1/2), the 
first mode boundary shifts to the left and the second mode boundary becomes 
the critical flutter boundary. Again, the results agree well with Cunninghaa's 


solution. 
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Figure 4. The first mode shape (real part) for panel 
corresponding to point A in Fig. 3 
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Figure 5 The second mods shipe (real part) for the 
center Jine of panel curresponding to 
point Uo in Fig. 3 
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The first and second flutter mode shapes (real part) corresponding 
to points A and B in Fig. 6 are plotted in Figs. 7 and 8, respectively. 


They are quite similar to those shown in Figs. 4 and 5&. 


Fig. 9 shows that, as the panel becomes square, the third and fourth 
mode boundaries emerge and the third mode boundary becomes the critical flutter 
boundary. The results are in good agreement with Cunningham's solution. 

The results for the first mode boundary obtained by Houbolt using the piston 


theory are also shown in the figure. 


The third and fourth flutter mode shapes (real part) corresponding to 
points A and B in Fig. 9 are shown in Figs. 10 and 11, respectively. The 
effect of the flow that tends to blow flat the front part of the panel is 


seen. 


Fig. 12 shows that, for a long panel with g/w = 2, the fifth mode flutter 
boundary becomes dominant in the low mass ratio range and the first mode 
boundary dominates the high mass ratio range. It also shows that Houbolt’s 
solution for first mode and this solution are very close. There are, however, 
some discrepancies between this solution and Cunningham's solution for the 
first mode boundary. This may be due to the fact that the stiffness coupling 


effect between assumed beam natural modes was uniformly neglected by Cunningham 


Such coupling effect can become significant as the chord-span ratio increases. 
It should be noted that, due to the limited number of elements and boxes used, 


the present results are not the most accurate that this approach can produce. 
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The third mode shape (reel part) for the 
center line of panel] corresponding to 
point A in Fig. 9 


Figure 13 
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The first and third flutter mode shapes (real part) corresponding to 
points A and B in Fig. 12 are shown in Figs. 13 and 14, respectively. The 
difference between a symmetrical in vacuo first mode and the first flutter 
mode with front part blown flat is clearly shown in ag. 13. 

Fig. 15 shows that, for a very long panel with ¢/w = 4, the tenth mode 
flutter boundary dominates the lower mass ratio range while the first mode 
boundary dominates the upper mass ratio range. The discrepancy between the 
present solution and Cunningham's solution may be attributed to the same reasons 


as explained for Fig. 6. Houbolt's solution for first mode boundary is also 


The second and tenth flutter mode shapes (real part) corresponding to 


points A and B in Fig. 15 are plotted in Figs. 16 and 17, respectively. 
Case 2 Clamped Panel with ¢/w = 2 at Various Mach Numbers 


A clausd panel with ¢/w = 2 and with one surface exposed to air stream 
with various Mach numbers was studied. One of the main purposes was to establish 
a curve for the thickness ratios required to prevent flutter of panel at various 


air speeds and at sea level. 


The first mode flutter boundaries for panel with */w = 2 and M = 1.05, 1.13 


1.4, 1.55 2.0, and 3,0 are shown in Figs. 18, 19, and 20, respectively. Corresponding 
tou each curve, a dashed parabola is shown, Lach parabola is plotted for the 


equatinn xy = © or (pe/onj{u 2/¥) = ©. The censtant C is dependent upon 
: 
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Figure 13. The first mode shape (real part) for panel 
corresponding to point A in Pig. 12. 


Figure 14 The fifth mode shape (real vart) for the 
center Line of panel correspondiny to 
potnt Bin Fig. 12. 
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the densities of the aluminum and the air at sea lever, the dimensions of the 
panel, the air speed, and the bending rigidity of the panel. The intersectina 
point between a pair of dashed and solid curves gives the thickness ratio 

h/2 required to prevent flutter of tne aluminum panel at sea level. For 

all Mach numbers considered here, the first mode flutter boundaries are 


critical for the dashed curves for the sea level altitude. 


The results are shown in Fig. 21. The present results are slightly on 


the unconservative side as compared to the results of Cunningham. 
Case 3 Pinned-Edge Panels with Various Chord-Snan Ratios at M = 1.3 


The panels with simply-supported edges and various chord length-span 
width ratios were studied. The flow speed considered was M = 1.3. The 
dominating flutter boundaries for panels with c/w = 0.75, 0.5, 1.0, and 2.0 are 
shown in Figs. 22, 24, 26 anc 28, respectively. They are all in good 


agreement with chose found by Cunningham (Reference 9). 


Fig. 22 shows that the first mode flutter boundary is the critical boundary 


for panels with Liw = 0.25. The mode shape (real part) is shown in Fig. 23. 


Fig, 24 shows that the second mode flutter boundary 15 the critical 
bounuary for panels with #/w = 0.5. The first and secona flutter mode 


shapes (real part) are shown in Fig. 25. 


Fig. 26 shows that the third mode flutter boundary is the critical 
boundary for panel with &/w = 1.0. The first mode boundary found hy 
Hedgeseth (Reference W) using the piston theory is also shown in tne fiaure. 


The corresponding first end third flutter mode shapes (real part) are stown 
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Fig. 28 shows that the fifth mode flutter boundary is the critical 
boundary for panels with £/w = 2.0. The results by Hedgepeth and 
Cunningham aré seen to be in good agreement with the present results. 
The corresponding first and fifth flutter mode shapes (real part) are 


shown in Fig. 29, 


Case 4 Clamped Square Panels with Various Tension Parameters and Mach Numbers 


One of the advantages of the finite element method is that the 
effect of in-plane forces can be included in a direct and accurate 
fashion. To demonstrate this, a square clamped panel was chosen and 
four different tension parameters F = 0.01, 0.7, 0.5, ana 1 were concidered. 


The flutter boundaries were obtained for various Mach numbers. 


Fig. 30 shows the dominating first mode flutter boundaries and 
the hyperbola for the square, aluminum panel at 25,000 feet above sea 


level and at M= 1.1. Four values of tension parameters were included. 


The mode shapes (real part) corresponding tc points A (F = 0.01) 
and point B (F = 0.1) of the flutter boundaries in Fig, 30 are plotted 
in Figs. 31 and 32. The mode shapes for the case F = 0.5 and 1.0 are 


similar to those shown in Figs. 31 and 32. They are thus not shown here. 


Fig. 33 shows the first and the third mode flutter boundaries 
for M = 1.2 and various tension parameters. For aluminum panels at 25,000 feet above 
sea level (dashed parabola), the third mode boundarics are the critical 
flutter boundaries. The first fluttcr mode shapes (real part) corresponding 


to points A and B in Fig. 33 are presented in Figs. 34 and 35, respectively. 
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The third flutter mode snapes (real part) for the center chord Vines 


corresponding to noints @, 0, E, and = are shown in Fig. 36. 


Fig. 37 shows the first mode flutter boundaries for M = 1.3 and the 
third mode flutter boundaries for M = 1.3, 1.32, 1.35 for various values 
of tension. The flutter boundaries are dominated by the third made. Since 
the tlutter mode shapes are similar to those shown in the previous figures, 


they are not presented here. 


Fig. 38 shows the Tiutter boundaries for M = 1.4, 1.45. 1.48 and 
various tension parameters. Tie third mode flutter boundaries shift te the 
right as the Mach number increases, It is interesting to see that the top 
portions of the third mode flutter boundarics begin to bend dowa to the left 


as the Mach number increases. 


Fig, 39 shows that for M = 1.6 and various values of tension, the 


third mode flutter boundaries continue to bend down rapidly with a slight 


increase in Mach miutaber. 


Figs. 40 and 41 show that for M = 1.52 and 1.54 and various v «Tues 
of tension, the third mode flutter buundaries bend down substantially to 


be small loops. They only dominaie the lower range of the mass ratio. 


Fig 42 shows that for M= 1.6, ali the third mode fiutter boundary ies 
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disappear and the first mode fiutter bounderies once again becone the critical 


boundaries. Such a phenomenon is also the case for M = 2.0 as shown in 


Fig. 43. : 
Also shown in Figs. 37-43 are the dashed paranolas for an aluminum pane: 


at 25,CQ0 feet above sea level. 
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Figure 40, Flutter boundaries for clamped panels of €/w=1.0 with M=1.52 
and various values of tension parameters F. 
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Figure 41 Flutter boundaries for clamped panels of £/w=1.0 with M=1.54 
and various values of tension parameters FPF, 
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Figure 42 Flutter boundaries for clamped panels of 2/w=1.0 with 
M=1.6 and various values of tension parameters F, 
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By collecting all the intersecting points between the dashed parabolas 
and the critical flutter boundaries in Figs. 30, 33, and 37-43, the results 
for thickness ratios required to prevent flutter of the panel are shown 


in Fig. 44 for various values of Mach numbers and tension. 


In Fig. 44, the critical flutte: houndaries were dominated by the 
third mode in the Mach region between approximately 1.2 and 1.5. For other 
lower and higher Mach regions, the first mode flutter boundaries -. ate. 
The sharp drops in the curves are due to the abrupt changes inci ~ ul 


flutter modes. 


The beneficial effect of introducing in-plane tensions to reduce the 


required panel thickness to avoid flutter is clearly demonstrated. Fig. 44 


should be of value to panel flutter analysts and designers. 
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SECTION V 
CONCLUDING REMARKS 


A basic finite element procedure for panel flutter analysis has been 
developed and performed with examples. The foliowing concluding remarks 


may be made. 


(1) The three-dimensional supersonic unsteady potential flow theory 
was employed. This theory allows the treatment of panels with finite aspect 
ratio. It is particularly advantageous at low supersonic range (1 < M< ¥2 
for panels with chievd-span ratio less than one, when the piston theory does 


not. give satisfactory results. 


(2) The finite clement method offers a set of elegant anda straightforward 
eigenvalue equations. It can be used directly to solve for flutter frequencies 
and mode shapes without requiring the natural vibration frequencies and mode 
Snapes before the cigensolution, and also without requiring the computation 


of mode shapes after the eigensolution. 


(3) Jf the modal method is used in panel flutter analysis, one cculd 
ase the natural frequencies and modes found by the finite element method. 
Such procedure is, however, basically different from the present finite 


element metho! 


(4) The furmuletions here (kgs. 19-24) are general. They are readily 
applicable to any plate finite element displacement model whose shape functions 


ere Kngon. 


ft mes Pn lil le a et a, 


(5) The present results agree well with Galerkin's modal solution 


by Cunningham. 


(6) The critical flutter boundaries for a square clamped panel 


change modes abruptly as the Mach number is varied. 


(7) Fig. 44 clearly demonstrates the beneficial effect of in-plane 


tensions. 


(8) The present development may provide a basic finite element procedure 
for panel flutter analysis. The present results may provide useful data to 


the flutter analysts and designers. 


(5) A jloyical next step appears to be the development of a method 
using triangular plate finite elements combined with triangular Mach boxes. 
Such development will be of great value to the flutter predictions of wing 


panels with arbitrary configurations. It is suggested that the triangular 


element developed by Bell (Re1erence 15) be used, 


he en a RN la ii ts atthe a 
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APPENDIX 
THE COMPUTER PROGRAM 


Ty. descripcion and listing of the computer program are provided in 
uwhis appendix. The flow chart is given on pages 69 and 70. The procedure 
for input is given on pages 62, 67, and 68. The use of the program is 
demonstrated by considering an example of a Square clamped panel with 


Ky = 0.5, M= 1.2, and F = 0.1. Both input and output data are provided. 


The stiffness, mass, and incremental stiffness matrices are based on 
the known explicit coefficients. These coefficients are read in as input 
data given on pages 63-66. The program uses the subroutine EISPACK for 


solving the genera! complex eigenvalue equations. 
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INPUT DATA 


un Control Card (2F5.2, 713, 4F7.4) ° 


Columns 1-5 
6-10 
11-13 
14-16 
17-19 


20-22 


23-25 


26-28 
29-31 
32-38 
39-45 
46-52 
53-59 


Panel chordwise length (XL) 

Panel spanwise width (YL) 

Number of elements in chordwise direction (NX) 
Half number of elements in spanwise direction (NY) 


Number of boxes in streamwise (or chordwise) direc- 
tion for one element (I6X) 


Humber of boxes in cross-streamwise (or spanwise) 
direction for one element (IBY) 


He) number of degrees of freedom for the panel 
(NXD 


Number of aodes in streamwise direction (NDX) 
Number of nodes in cross-streamwise direction (NDY) 
Structural damping coefficient (SG) 

Mach number (MACH) 

Starting value for the mass ratio range (MUO) 


Increment for the values of mass ratios (XINCR) 


2. Element matrix cards (814, 212, 214, 216) 
Read in the data given on pages 63-66. 
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DATA INFORMATION FOR QQ(I.J), QM(I,J) and Qu(I,d), continued ... 
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Element Cards with the boundary condition code at each node (2213). 


Columns 


La3 
4-6 


7-9 


10~12 


13-15 


16-18 


19-2] 


34-36 


37-39 


40-42 


The sequence number for the element 


The starting matrix row position for the 
degrees of freedom at nodal point I 


The starting matrix row position for the 
degrees of freedom at nudal peint J 


The Starting matrix row position for the 
degrees of freedom at nedal point kK 


The starting nmatrix row position for the 
degrees of freedom at nodal point L 


Boundary condition code for deflection 
Z at nodal point I 


Boundary condition code for Z . at nodal 


point I : 

Boundary condition code for z y at nodal 
point I , 

Boundary condition code for 2 sd at nodal 
point I ay, 


Boundary condition code for Z at nodal 
point J 


Boundary condition code ror Z | at nodal 
eee Pe 
poink uv 


Boundary condition cede tor Z y at nodal 


point J , 

Boundary condition code for Z at nodal 
iy xy 

point J “ 


Boundary condition coce for Z at nodal 
point K 


Boundary condition code for 2 i at nodal 
point kK 1 


Goundary condition coce for Z | at nodal 
point ¥ a 


Boundary candition cade for Z a, at nodal 
point \ ve 


6) 


(1D(1)) 


52-54 Boundary condition code for Z at nodal 
point L (w(13)) 


55-57 Boundary condition code for Z . at nodal 
point..L : (w(14)) 


58-60 Boundary condition code for Z at nodal 
point L 2 (w(15)) 


61-63 Boundary condition code for Z at nodal 
point L oXy (w(16)) 


The 'ID' array at I, J, K and L is assigned for each element for 
a specific mesh of the elements of the panel. The permutation for 1, J, 
K and L is clackwise and I is corresponding to the origin of each local 
element as shown in Fig. 2. If ID(I) = 0, all the degrees of freedom at 
node I are zero, Same situation is apnlied for other nodes for the 
element. 

The ‘w' array defines the boundary condition codes tor the degrees of 
freedom at each node. If w(t) = 0, the associated degree of freedom is 
zero and if w(I) # 0, the associated degree ef freedom is taken into 


account. 
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